SUBROUTINE WRTPCDF
! (C) Copyright 2000- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!**** *WRTPCDF*  - Writing prognostic variables of the one-column surface model

!     Purpose.
!     --------
!     Write out prognostic variables in NetCDF format

!**   Interface.
!     ----------
!        *CALL* *WRTPCDF

!        Explicit arguments :
!        --------------------


!        Implicit arguments :
!        --------------------

!     Method.
!     -------
!        See documentation

!     Externals.
!     ----------
!        None

!     Reference.
!     ----------
!        ECMWF Research Department documentation 
!        of the one-column surface model

!     Author.
!     -------
!        Bart vd Hurk *KNMI*

!     Modifications.
!     --------------
!        Original : 2000-07-24
!        Y. Takaya (ECMWF) Add ocean mixed layer model output 2008-10-07

!     ------------------------------------------------------------------


USE YOMLUN1S , ONLY : NULOUT   ,NPOSGG   ,NPOSOCP,  RMISS
USE YOMCT01S , ONLY : NSTART   ,NFRPOS
USE YOMDYN1S , ONLY : NSTEP    ,TSTEP
USE YOETHF   , ONLY : RHOH2O
USE YOMLOG1S , ONLY : LDBGS1   ,CFFORC   ,NACCUR   ,NDIMCDF  ,&
            &LWRGG,   LWROCP,LNCSNC
USE YOMGP1S0 , ONLY : TSLNU0   ,QLINU0   ,TILNU0   ,FSNNU0   ,&
            &TSNNU0   ,ASNNU0   ,RSNNU0   ,TRENU0   ,WRENU0  ,& 
            &TLICENU0,TLMNWNU0,TLWMLNU0,TLBOTNU0,TLSFNU0,& ! ENDUTRA
            &HLICENU0,HLMLNU0 ,WSNNU0,&                           ! ENDUTRA 
            &UONU0    ,VONU0    ,TONU0    ,SONU0           ! KPP
USE YOMDPHY  , ONLY : NLAT     ,NLON     ,NPOI     ,    NPOIOFF    ,&
            &NLALO   ,NCSS     ,&
            &NCOM,NCSNEC,NPOIP                                          !KPP
USE YOESOIL1S, ONLY : RDAW
USE YOMGC1S  , ONLY : LMASK
USE YOMRIP   , ONLY : RSTATI
USE PARKIND1 ,ONLY : JPIM     ,JPRB     ,JPRD,  JPRM
USE YOMHOOK  ,ONLY : LHOOK    ,DR_HOOK, JPHOOK 
USE YOEPHY   , ONLY : LEFLAKE , LESNML
USE YOMGDI1S , ONLY : D1STITK2
USE MPL_MODULE
IMPLICIT NONE

!* Netcdf Interface
INTEGER ISTART1(1),ISTART3(3),ISTART4(4)
INTEGER ICOUNT1(1),ICOUNT3(3),ICOUNT4(4)
REAL(KIND=JPRM),ALLOCATABLE :: ZOUTPUTS(:,:)
REAL(KIND=JPRD),ALLOCATABLE :: ZOUTPUT(:,:)
REAL(KIND=JPRD) ZTIME
INTEGER NSTPP,NPOS,NVARID,IERR

!* Local variables
INTEGER(KIND=JPIM), PARAMETER :: JPNCDF=2
INTEGER(KIND=JPIM) :: IPOS(JPNCDF)
LOGICAL LPOS(JPNCDF)
INTEGER(KIND=JPIM) :: IX1,IX2,IY1,IY2,NLEN,IST,IEN,JSLEV,JK,JL,&
     &    IYYMD,IHM,J,ISTP,IENP
REAL(KIND=JPRB) :: RTT,RLMLT,ZIHCAP,ZEPSILON
REAL(KIND=JPRD) :: ZJUL

REAL(KIND=JPRB),ALLOCATABLE :: ZQLINU0(:,:),ZSSN(:),ZRSN(:),ZTSN(:),ZWSN(:),ZTMP(:),ZTMP1(:)
REAL(KIND=JPRB),ALLOCATABLE :: ZBUF(:)
INTEGER(KIND=JPIM) :: MYPROC, NPROC
INTEGER(KIND=JPIM) :: NMXNV
REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

#include "netcdf.inc"

IF (LHOOK) CALL DR_HOOK('WRTPCDF',0,ZHOOK_HANDLE)

MYPROC = MPL_MYRANK()
NPROC  = MPL_NPROC()

ISTP   = NPOIOFF(1)+1
IENP   = NPOIOFF(NPROC+1)

WRITE(NULOUT,*) 'WRTPCDF'


ALLOCATE (ZQLINU0(NPOI,NCSS))  !RESCALED SOIL WATER
ALLOCATE (ZSSN(NPOI))
ALLOCATE (ZRSN(NPOI))
ALLOCATE (ZTSN(NPOI))
ALLOCATE (ZWSN(NPOI))
ALLOCATE (ZTMP(NPOI))
ALLOCATE (ZTMP1(NPOI))

ZEPSILON=10._JPRB*EPSILON(ZEPSILON)

IF(NDIMCDF == 2)THEN
  IX1=1
  IX2=NLON
  IY1=1
  IY2=NLAT
  NLEN=(IX2-IX1+1)*(IY2-IY1+1)
  IST=IX1+(IY1-1)*NLON
  IEN=IST+NLEN-1
ELSE
  IST=1
  IEN=NLALO
  NLEN=IEN-IST+1
ENDIF

NMXNV=MAX(NCSS,(MAX(NCSNEC,NCOM+1)))

ALLOCATE (ZOUTPUT(NLEN,NMXNV))        !INFLATED OUTPUT ARRAY
ALLOCATE (ZOUTPUTS(NLEN,NMXNV))           !INFLATED OUTPUT ARRAY

ALLOCATE(ZBUF(NLEN))

ZOUTPUT=RMISS
ZBUF=RMISS
!     ------------------------------------------------------------------

!*       1.   SCALING SOME VARIABLES.
!             -----------------------

!     NOTE: UNITS OF FSNNU,WRENU IN THE SURF1D MODEL:            KG/M2
!           UNITS OF SNOW AND INTERCEP RESER WATER IN THE
!                                     PHYSICS:                   KG/M2
!           UNITS OF SNOW AND INTERCEP RESER WATER IN THE
!                                     INPUT AND OUTPUT FILES:    KG/M2


!     NOTE: UNITS OF QLINU IN THE SURF1D MODEL:                  M3 M-3
!           UNITS OF SOIL WATER IN THE INPUT FILE:               M3 M-3
! UNITS OF SOIL WATER IN THE INPUT FILE IF RESTART:              KG/M2/LAYER
!           UNITS OF QLINU IN THE OUTPUT:                        KG/M2/LAYER


DO JK=1,NCSS
  ZQLINU0(1:NPOI,JK)=QLINU0(1:NPOI,JK)*RDAW(JK)*RHOH2O
ENDDO

RTT=273.16_JPRB
RLMLT=2.8345E+6_JPRB-2.5008E+6_JPRB
ZIHCAP=2.05E6_JPRB/920._JPRB
! Total snow variables
ZSSN(:) = SUM(FSNNU0(:,:),DIM=2)
ZTMP(:) = SUM(FSNNU0(:,:)/RSNNU0(:,:),DIM=2) ! SNOW DEPTH 
ZTMP1(:) = SUM(ZIHCAP*FSNNU0(:,:)*(TSNNU0(:,:)-RTT) - &    ! INTERNAL ENERGY
              RLMLT*(FSNNU0(:,:)-WSNNU0(:,:)),DIM=2)
WHERE (ZSSN > ZEPSILON .AND. ZTMP > ZEPSILON )
  ZRSN = ZSSN / ZTMP
  ZTSN = MIN(RTT, RTT + (ZTMP1+RLMLT*ZSSN)/(ZIHCAP*ZSSN))
!!ZTSN = MIN(RTT, TSNNU0(:,1))
  ZWSN = MAX(0._JPRB,ZSSN+ZTMP1/RLMLT)
ELSEWHERE
  ZRSN = 100._JPRB
  ZTSN = RTT
  ZWSN = 0._JPRB
ENDWHERE


DEALLOCATE(ZTMP1)
DEALLOCATE(ZTMP)


!*       2.    WRITE TIME STEP.
!              ----------------

IF (LDBGS1) THEN
  JL=1
  WRITE(NULOUT,*) 'TIME STEP'
  WRITE(NULOUT,'(I6)') NSTEP

  WRITE(NULOUT,*) 'PROGNOSTIC VARIABLES'
  WRITE(NULOUT,*) 'GRID BOX JL = ',JL

!     -----------------------------------------------------------------


!*       3.    WRITE SOIL VARIABLES.
!              ---------------------

  WRITE(NULOUT,*) ' SOIL TEMPERATURE - SOIL MOISTURE - ICE T'
  DO JSLEV=1,NCSS
    WRITE(NULOUT,'(3(2X,E12.6))') TSLNU0(JL,JSLEV),&
     &QLINU0(JL,JSLEV),&
     &TILNU0(JL,JSLEV)
  ENDDO


!*       4.    WRITE SKIN VARIABLES.
!              ---------------------

  WRITE(NULOUT,*) ' SKIN TEMP. - SKIN. RES. CONT. '
  WRITE(NULOUT,'(2(2X,E12.6))') TRENU0(JL) , WRENU0(JL)


!*       5.    WRITE SNOW VARIABLES
!              -----------------

  WRITE(NULOUT,*) 'SNOW_DEPTH   SNOW_T  SNOW_ALBEDO  SNOW_DENSITY'
  WRITE(NULOUT,'(E12.6)') FSNNU0(JL,1),TSNNU0(JL,1),ASNNU0(JL),RSNNU0(JL,1)
ENDIF

!        6. WRITE TO NETCDF FILES


IPOS(1)=NPOSGG
IPOS(2)=NPOSOCP
LPOS(1)=LWRGG
LPOS(2)=LWROCP

!     time level
CALL DATTIM(zjul,IYyMD,IHM)
NSTPP=NSTEP
ZTIME=REAL(RSTATI,KIND=JPRD)-REAL(0.5_JPRD*TSTEP,KIND=JPRD)
ISTART1(1)=(NSTEP-NSTART)/NFRPOS+1
ICOUNT1(1)=1
IF( MYPROC == 1 ) THEN
  DO J=1,JPNCDF
    IF(LPOS(J))THEN
      NPOS = IPOS(J)
      NVARID = NCVID(NPOS,'time',IERR)
      CALL NCVPT (NPOS,NVARID,ISTART1,ICOUNT1,ZTIME,IERR)
      NVARID = NCVID(NPOS,'timestp',IERR)
      CALL NCVPT (NPOS,NVARID,ISTART1,ICOUNT1,NSTPP,IERR)
    ENDIF
  ENDDO
ENDIF

IF(NDIMCDF == 2)THEN
  ISTART3(1)=IX1
  ISTART3(2)=IY1
  ISTART3(3)=(NSTEP-NSTART)/NFRPOS+1
  ICOUNT3(1)=IX2-IX1+1
  ICOUNT3(2)=IY2-IY1+1
  ICOUNT3(3)=1
ELSE
  ISTART3(1)=IST
  ISTART3(2)=(NSTEP-NSTART)/NFRPOS+1
  ICOUNT3(1)=NLEN
  ICOUNT3(2)=1
ENDIF
      
IF(LWRGG)THEN

!*    THE PROGNOSTIC OUTPUT FILE
      
  NPOS=NPOSGG
      
!     skin temperature
  CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=TRENU0(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:AvgSurfT")
  IF( MYPROC == 1 ) THEN
    NVARID = NCVID(NPOS,'AvgSurfT',IERR)
    ZOUTPUT(:,1)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
    IF(NACCUR == 1)THEN
      ZOUTPUTS(:,1)=ZOUTPUT(:,1)
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUTS(:,1),IERR)
    ELSE
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT(:,1),IERR)
    ENDIF
  ENDIF
        
!     skin layer content
  CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=WRENU0(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:CanopInt")
  IF( MYPROC == 1 ) THEN
    NVARID = NCVID(NPOS,'CanopInt',IERR)
    ZOUTPUT(:,1)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
    IF(NACCUR == 1)THEN
      ZOUTPUTS(:,1)=ZOUTPUT(:,1)
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUTS(:,1),IERR)
    ELSE
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT(:,1),IERR)
    ENDIF
  ENDIF

!     snow water equivalent
  CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZSSN(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:SWE")
  IF( MYPROC == 1 ) THEN
    NVARID = NCVID(NPOS,'SWE',IERR)
    ZOUTPUT(:,1)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
    IF(NACCUR == 1)THEN
      ZOUTPUTS(:,1)=ZOUTPUT(:,1)
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUTS(:,1),IERR)
    ELSE
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT(:,1),IERR)
    ENDIF
  ENDIF

  DEALLOCATE(ZSSN)

!     snow temp
  CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZTSN(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:SnowT")
  IF( MYPROC == 1 ) THEN
    NVARID = NCVID(NPOS,'SnowT',IERR)
    ZOUTPUT(:,1)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
    IF(NACCUR == 1)THEN
      ZOUTPUTS(:,1)=ZOUTPUT(:,1)
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUTS(:,1),IERR)
    ELSE
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT(:,1),IERR)
    ENDIF
  ENDIF
 
  DEALLOCATE(ZTSN)
 
!     snow temp
  CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZWSN(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:slw")
  IF( MYPROC == 1 ) THEN
    NVARID = NCVID(NPOS,'slw',IERR)
    ZOUTPUT(:,1)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
    IF(NACCUR == 1)THEN
      ZOUTPUTS(:,1)=ZOUTPUT(:,1)
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUTS(:,1),IERR)
    ELSE
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT(:,1),IERR)
    ENDIF
  ENDIF

  DEALLOCATE(ZWSN)

!     snow density
  CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZRSN(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:snowdens")
  IF( MYPROC == 1 ) THEN
    NVARID = NCVID(NPOS,'snowdens',IERR)
    ZOUTPUT(:,1)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
    IF(NACCUR == 1)THEN
      ZOUTPUTS(:,1)=ZOUTPUT(:,1)
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUTS(:,1),IERR)
    ELSE
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT(:,1),IERR)
    ENDIF
  ENDIF

  DEALLOCATE(ZRSN)

!     snow albedo
  CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ASNNU0(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:SAlbedo")
  IF( MYPROC == 1) THEN
    NVARID = NCVID(NPOS,'SAlbedo',IERR)
    ZOUTPUT(:,1)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
    IF(NACCUR == 1)THEN
      ZOUTPUTS(:,1)=ZOUTPUT(:,1)
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUTS(:,1),IERR)
    ELSE
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT(:,1),IERR)
    ENDIF
  ENDIF
  
!* < ENDUTRA 
  IF (LEFLAKE) THEN
    ! lake ice temperature
    CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=TLICENU0(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:TLICE")
    IF( MYPROC == 1 ) THEN
      NVARID = NCVID(NPOS,'TLICE',IERR)
      ZOUTPUT(:,1)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
      IF(NACCUR == 1)THEN
        ZOUTPUTS(:,1)=ZOUTPUT(:,1)
        CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUTS(:,1),IERR)
      ELSE
        CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT(:,1),IERR)
      ENDIF
    ENDIF
  
    !     lake mean water temperature
    CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=TLMNWNU0(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:TLMNW")
    IF( MYPROC == 1 ) THEN
      NVARID = NCVID(NPOS,'TLMNW',IERR)
      ZOUTPUT(:,1)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
      IF(NACCUR == 1)THEN
        ZOUTPUTS(:,1)=ZOUTPUT(:,1)
        CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUTS(:,1),IERR)
      ELSE
        CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT(:,1),IERR)
      ENDIF
    ENDIF
  
    !     lake mixed layer temperature
    CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=TLWMLNU0(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:TLWML")
    IF( MYPROC == 1 ) THEN
      NVARID = NCVID(NPOS,'TLWML',IERR)
      ZOUTPUT(:,1)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
      IF(NACCUR == 1)THEN
        ZOUTPUTS(:,1)=ZOUTPUT(:,1)
        CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUTS(:,1),IERR)
      ELSE
        CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT(:,1),IERR)
      ENDIF
    ENDIF
  
    !     lake bottom temperature 
    CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=TLBOTNU0(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:TLBOT")
    IF( MYPROC == 1 ) THEN
      NVARID = NCVID(NPOS,'TLBOT',IERR)
      ZOUTPUT(:,1)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
      IF(NACCUR == 1)THEN
        ZOUTPUTS(:,1)=ZOUTPUT(:,1)
        CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUTS(:,1),IERR)
      ELSE
        CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT(:,1),IERR)
      ENDIF
    ENDIF
  
    !     lake temperature shape factor
    CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=TLSFNU0(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:TLSF")
    IF( MYPROC == 1 ) THEN
      NVARID = NCVID(NPOS,'TLSF',IERR)
      ZOUTPUT(:,1)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
      IF(NACCUR == 1)THEN
        ZOUTPUTS(:,1)=ZOUTPUT(:,1)
        CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUTS(:,1),IERR)
      ELSE
        CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT(:,1),IERR)
      ENDIF
    ENDIF
  
    !     lake ice thickness
    CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=HLICENU0(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:HLICE")
    IF( MYPROC == 1 ) THEN
      NVARID = NCVID(NPOS,'HLICE',IERR)
      ZOUTPUT(:,1)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
      IF(NACCUR == 1)THEN
        ZOUTPUTS(:,1)=ZOUTPUT(:,1)
        CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUTS(:,1),IERR)
      ELSE
        CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT(:,1),IERR)
      ENDIF
    ENDIF
  
    !     lake mixed layer thickness
    CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=HLMLNU0(:),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:HLML")
    IF( MYPROC == 1 ) THEN
      NVARID = NCVID(NPOS,'HLML',IERR)
      ZOUTPUT(:,1)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
      IF(NACCUR == 1)THEN
        ZOUTPUTS(:,1)=ZOUTPUT(:,1)
        CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUTS(:,1),IERR)
      ELSE
        CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT(:,1),IERR)
      ENDIF
    ENDIF
  ENDIF

!* ENDUTRA > 

  IF(NDIMCDF == 2)THEN
    ISTART4(1)=IX1
    ISTART4(2)=IY1
    ISTART4(3)=1
    ISTART4(4)=(NSTEP-NSTART)/NFRPOS+1
    ICOUNT4(1)=IX2-IX1+1
    ICOUNT4(2)=IY2-IY1+1
    ICOUNT4(3)=NCSS
    ICOUNT4(4)=1
  ELSE
    ISTART4(1)=IST
    ISTART4(2)=1
    ISTART4(3)=(NSTEP-NSTART)/NFRPOS+1
    ICOUNT4(1)=NLEN
    ICOUNT4(2)=NCSS
    ICOUNT4(3)=1
  ENDIF
!     soil temperature
  IF( MYPROC == 1 ) NVARID = NCVID(NPOS,'SoilTemp',IERR)
  DO JK=1,NCSS
    CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=TSLNU0(:,JK),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:SoilTemp")
    IF( MYPROC == 1 ) THEN
      ZOUTPUT(:,JK)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
    ENDIF
  ENDDO
  
  IF( MYPROC == 1 ) THEN
    IF(NACCUR == 1)THEN
      ZOUTPUTS(:,1:NCSS)=ZOUTPUT(:,1:NCSS)
      CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUTS(:,1:NCSS),IERR)
    ELSE
      CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUT(:,1:NCSS),IERR)
    ENDIF
  ENDIF

!     soil moisture
  IF( MYPROC == 1 ) NVARID = NCVID(NPOS,'SoilMoist',IERR)
  DO JK=1,NCSS
    CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZQLINU0(:,JK),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:SoilMoist")
    IF( MYPROC == 1 ) THEN
      ZOUTPUT(:,JK)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
    ENDIF
  ENDDO

  DEALLOCATE(ZQLINU0)

  IF( MYPROC == 1 ) THEN
    IF(NACCUR == 1)THEN
      ZOUTPUTS(:,1:NCSS)=ZOUTPUT(:,1:NCSS)
      CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUTS(:,1:NCSS),IERR)
    ELSE
      CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUT(:,1:NCSS),IERR)
    ENDIF
  ENDIF 

!     sea ice temperature
  IF( MYPROC == 1 ) NVARID = NCVID(NPOS,'icetemp',IERR)
  DO JK=1,NCSS
    CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=TILNU0(:,JK),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:icetemp") 
    IF( MYPROC == 1 ) THEN
      ZOUTPUT(:,JK)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
    ENDIF
  ENDDO

  IF( MYPROC == 1 ) THEN
    IF(NACCUR == 1)THEN
      ZOUTPUTS(:,1:NCSS)=ZOUTPUT(:,1:NCSS)
      CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUTS(:,1:NCSS),IERR)
    ELSE
      CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUT(:,1:NCSS),IERR)
    ENDIF
  ENDIF
  
  IF (LESNML) THEN
    IF(NDIMCDF == 2)THEN
      ICOUNT4(3)=NCSNEC
    ELSE
      ICOUNT4(2)=NCSNEC
    ENDIF
    
    !     SWE
    IF( MYPROC == 1 ) NVARID = NCVID(NPOS,'SWEML',IERR)
    DO JK=1,NCSNEC
      CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=FSNNU0(:,JK),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:SWEML") 
      IF( MYPROC == 1 ) THEN
        ZOUTPUT(:,JK)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
      ENDIF
    ENDDO

    IF( MYPROC == 1 ) THEN
      IF(NACCUR == 1)THEN
        ZOUTPUTS(:,1:NCSNEC)=ZOUTPUT(:,1:NCSNEC)
        CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUTS(:,1:NCSNEC),IERR)
      ELSE
        CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUT(:,1:NCSNEC),IERR)
      ENDIF
    ENDIF
    !     SwnoT
    IF( MYPROC == 1 ) NVARID = NCVID(NPOS,'SnowTML',IERR)
    DO JK=1,NCSNEC
      CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=TSNNU0(:,JK),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:SnowTML") 
      IF( MYPROC == 1 ) THEN
        ZOUTPUT(:,JK)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
      ENDIF
    ENDDO

    IF( MYPROC == 1 )THEN
      IF(NACCUR == 1)THEN
        ZOUTPUTS(:,1:NCSNEC)=ZOUTPUT(:,1:NCSNEC)
        CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUTS(:,1:NCSNEC),IERR)
      ELSE
        CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUT(:,1:NCSNEC),IERR)
      ENDIF
    ENDIF

    !     Snowdensity 
    IF( MYPROC == 1 ) NVARID = NCVID(NPOS,'snowdensML',IERR)
    DO JK=1,NCSNEC
      CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=RSNNU0(:,JK),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:snowdensML") 
      IF( MYPROC == 1 ) THEN
        ZOUTPUT(:,JK)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
      ENDIF
    ENDDO

    IF( MYPROC == 1 ) THEN
      IF(NACCUR == 1)THEN
        ZOUTPUTS(:,1:NCSNEC)=ZOUTPUT(:,1:NCSNEC)
        CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUTS(:,1:NCSNEC),IERR)
      ELSE
        CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUT(:,1:NCSNEC),IERR)
      ENDIF
    ENDIF

    !     Snow liquid water 
    IF( MYPROC == 1 ) NVARID = NCVID(NPOS,'slwML',IERR)
    DO JK=1,NCSNEC
      CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=WSNNU0(:,JK),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:slwML")
      IF( MYPROC == 1 ) THEN 
        ZOUTPUT(:,JK)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
      ENDIF
    ENDDO

    IF( MYPROC == 1 ) THEN
      IF(NACCUR == 1)THEN
        ZOUTPUTS(:,1:NCSNEC)=ZOUTPUT(:,1:NCSNEC)
        CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUTS(:,1:NCSNEC),IERR)
      ELSE
        CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUT(:,1:NCSNEC),IERR)
      ENDIF
    ENDIF
  ENDIF 
  
! synchronise
  1000 continue
  IF( MYPROC == 1 ) THEN
    IF (LNCSNC) CALL NCSNC(NPOS,IERR)
  ENDIF
ENDIF

!--Ocean mixed layer
IF(LWROCP)THEN

!*    THE PROGNOSTIC OUTPUT FILE
      
  NPOS=NPOSOCP

  IF(NDIMCDF == 2)THEN
    ISTART3(1)=IX1
    ISTART3(2)=IY1
    ISTART3(3)=(NSTEP-NSTART)/NFRPOS+1
    ICOUNT3(1)=IX2-IX1+1
    ICOUNT3(2)=IY2-IY1+1
    ICOUNT3(3)=1
  ELSE
    ISTART3(1)=IST
    ISTART3(2)=(NSTEP-NSTART)/NFRPOS+1
    ICOUNT3(1)=NLEN
    ICOUNT3(2)=1
  ENDIF
  
!     sea surface temperature
  CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=TONU0(:,1),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:OceanSST")
  IF( MYPROC == 1 ) THEN
    NVARID = NCVID(NPOS,'OceanSST',IERR)
    ZOUTPUT(:,1)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
    IF(NACCUR == 1)THEN
      ZOUTPUTS(:,1)=ZOUTPUT(:,1)
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUTS(:,1),IERR)
    ELSE
      CALL NCVPT (NPOS,NVARID,ISTART3,ICOUNT3,ZOUTPUT(:,1),IERR)
    ENDIF
  ENDIF
    
  IF(.TRUE.) THEN ! output ocean 3D fields?

    IF(NDIMCDF == 2)THEN
      ISTART4(1)=IX1
      ISTART4(2)=IY1
      ISTART4(3)=1
      ISTART4(4)=(NSTEP-NSTART)/NFRPOS+1
      ICOUNT4(1)=IX2-IX1+1
      ICOUNT4(2)=IY2-IY1+1
      ICOUNT4(3)=NCOM+1
      ICOUNT4(4)=1
    ELSE
      ISTART4(1)=IST
      ISTART4(2)=1
      ISTART4(3)=(NSTEP-NSTART)/NFRPOS+1
      ICOUNT4(1)=NLEN
      ICOUNT4(2)=NCOM+1
      ICOUNT4(3)=1
    ENDIF

!     Ocean Temperature
    IF( MYPROC == 1 ) NVARID = NCVID(NPOS,'OceanT',IERR)
    DO JK=1,NCOM+1
      CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=TONU0(:,JK),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:OceanT")
      IF( MYPROC == 1 ) THEN
        ZOUTPUT(:,JK)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
      ENDIF
    ENDDO
    
    IF( MYPROC == 1 ) THEN
      IF(NACCUR == 1)THEN
        ZOUTPUTS(:,1:NCOM+1)=ZOUTPUT(:,1:NCOM+1)
        CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUTS(:,1:NCOM+1),IERR)
      ELSE
        CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUT(:,1:NCOM+1),IERR)
      ENDIF
    ENDIF

!     Ocean Salinity
    IF( MYPROC == 1 ) NVARID = NCVID(NPOS,'OceanS',IERR)
    DO JK=1,NCOM+1
      CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=SONU0(:,JK),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:OceanS")
      IF( MYPROC == 1 ) THEN
        ZOUTPUT(:,JK)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
      ENDIF
    ENDDO

    IF( MYPROC == 1 ) THEN
      IF(NACCUR == 1)THEN
        ZOUTPUTS(:,1:NCOM+1)=ZOUTPUT(:,1:NCOM+1)
        CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUTS(:,1:NCOM+1),IERR)
      ELSE
        CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUT(:,1:NCOM+1),IERR)
      ENDIF
    ENDIF

!     Ocean U velocity
    IF( MYPROC == 1 ) NVARID = NCVID(NPOS,'OceanU',IERR)
    DO JK=1,NCOM+1
      CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=UONU0(:,JK),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:OceanU")
      IF( MYPROC == 1 ) THEN
        ZOUTPUT(:,JK)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
      ENDIF
    ENDDO

    IF( MYPROC == 1 ) THEN
      IF(NACCUR == 1)THEN
        ZOUTPUTS(:,1:NCOM+1)=ZOUTPUT(:,1:NCOM+1)
        CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUTS(:,1:NCOM+1),IERR)
      ELSE
        CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUT(:,1:NCOM+1),IERR)
      ENDIF
    ENDIF

!     Ocean V velocity
    IF( MYPROC == 1 ) NVARID = NCVID(NPOS,'OceanV',IERR)
    DO JK=1,NCOM+1
      CALL MPL_GATHERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=VONU0(:,JK),KRECVCOUNTS=NPOIP(:),CDSTRING="WRTPCDF:OceanV")
      IF( MYPROC == 1 ) THEN
        ZOUTPUT(:,JK)=UNPACK(ZBUF(:),LMASK(ISTP:IENP),RMISS)
      ENDIF
    ENDDO

    IF( MYPROC == 1 ) THEN
      IF(NACCUR == 1)THEN
        ZOUTPUTS(:,1:NCOM+1)=ZOUTPUT(:,1:NCOM+1)
        CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUTS(:,1:NCOM+1),IERR)
      ELSE
        CALL NCVPT (NPOS,NVARID,ISTART4,ICOUNT4,ZOUTPUT(:,1:NCOM+1),IERR)
      ENDIF
    ENDIF

  ENDIF

!   CALL NCSNC(NPOS,IERR)
ENDIF

DEALLOCATE (ZOUTPUT)
DEALLOCATE (ZOUTPUTS)
DEALLOCATE (ZBUF)

IF (LHOOK) CALL DR_HOOK('WRTPCDF',1,ZHOOK_HANDLE)
      
RETURN

END SUBROUTINE WRTPCDF
